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In this paper, we present an efficient algorithm for the certification of numeric solu- 
tions to eigenproblems. The algorithm relies on a mixture of ball arithmetic, a suitable 
Newton iteration, and clustering of eigenvalues that are close. 
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1. INTRODUCTION 


Let IF be the set of floating point numbers for a fixed precision and a fixed exponent range. 
We will denote F? = {x € F: x 2 0). Consider an n x n matrix M € FF[i|l *" with complex 
floating entries. The numeric eigenproblem associated to M is to compute a transformation 
matrix T € IF[i|^ *" and a diagonal matrix D c IF[i]|" *" such that 


D x T MT. (1) 


The entries of D are the approximate eigenvalues and the columns of T' are the approximate 
eigenvectors of M. In addition, we might require that T is normalized. For instance, each 
of the columns might have unit norm. Alternatively, the norm of the i-th column may be 
required to be the same as the norm of the i-th row of T'^!, for each i. There are several 
well-known algorithms for solving the numeric eigenproblem [4]. 

Unfortunately, (1) is only an approximate equality. It is sometimes important to have 
rigourous bounds for the distance between the approximate eigenvalues and/or eigenvectors 
and the genuine ones. More precisely, we may ask for a diagonal matrix D, € (IF7)"*" and 
a matrix T; € (IF7)"*" of radii such that there exists a matrix T” € C"** for which 


D = ar 
is diagonal and 
[Di;— Dil 


[d 


(D): 


< 
< (T) 


n 
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for all i, j. This task will be called the certification problem of the numeric solution (D, T) 
to the eigenproblem for M. 

It will be convenient to rely on ball arithmetic |6, 8], which is a systematic technique 
for bound computations. When computing with complex numbers, ball arithmetic is more 
accurate than more classical interval arithmetic [11, 1, 12, 7, 9, 15], especially in multiple 
precision contexts. We will write IB — B(IF[i], F?) for the set of balls z = B(zc, 2,) = {z € €: 
|z — zel < zr} with centers ze in F[i] and radii z, in F”. In a similar way, we may consider 
matricial balls M — B(M., M,) € B(IF[i]" *^, (IF2)"*"): given a center matrix Me € F[i]”*” 
and a radius matrix M, € (F7)"*", we have 


M = B(M., M,) = {M EC": Vi, j, Mc) — Mij| < (Mr)uj]- 


Alternatively, we may regard B(M., M,) as the set of matrices in IB"*" with ball coeffi- 
cients: 


B(M., M; )i, j = B((Mc)i,j, (Mr)i,5): 


Standard arithmetic operations on balls are carried out in a reliable way. For instance, if u, 
v € B, then the computation of the product w = u v using ball arithmetic has the property 
that uv € w for any u € u and v € v. Given a ball z € B, it will finally be convenient to 
write |z| € F? and [z] € IF? for certified lower and upper bounds of |z| in F>. 


In the language of ball arithmetic, it is natural to allow for small errors in the input 
and replace the numeric input M € FF[i]|"*" by a ball input B(M., M,) € B”*”. Then we 
may still compute a numeric solution 


De ~ T MT, (2) 


for the eigenproblem associated to the center Me. Assume that the matrices in B(Me, My) 
are all diagonalizable. The generalized certification problem now consists of the compu- 
tation of a diagonal matrix D, € (IF7)"*" and a matrix T; € F|i]”*” such that, for every 
M € B(M., Mr), there exist D € B(D., Dr) and T € B(T;, Tj) with 


D = TMT. 


In absence of multiple eigenvalues, known algorithms for solving this problem such as |17, 
14] proceed by the individual certification of each eigenvector, which results in an O(n*) 
running time. 

Extensions to a cluster of eigenvalues and the corresponding eigenvectors have been 
considered in [3, 16], with similar O(n*) complexity bounds. Fixed points theorem based 
on interval arithmetic are used to prove the existence of a matrix with a given Jordan block 
in the matrix interval domain. Such an approach has been exploited for the analysis of 
multiple roots in [5, 13]. A test that provides an enclosing of all the eigenvalues has been 
proposed in [10]. Its certification relies on interval and ball arithmetics. The complexity 
of the test is in O(n?) but no iteration converging to the solution of the eigenproblem is 
described. 


In this paper, we present a new algorithm of time complexity O(n?) for certifying and 
enclosing clusters of eigenvectors and eigenvalues in a single step. We also provide an 
iterative procedure that converges geometrically to clusters of solutions. This convergence 
is quadratic in the case of single eigenvalues. Our algorithm extends a previous algorithm 
from [6] to the case of multiple eigenvalues. This yields an efficient test for approximate 
eigenvalues in the sense of the a-theory [2]. 
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We recall that it is very unlikely that the numeric matrix Me € IF[i|"*" with complex 
floating point coefficients has multiple eigenvalues. Indeed, small perturbations of matrices 
with multiple eigenvalues, as induced by rounding errors, generically only have simple 
eigenvalues. Consequently, we may assume without loss of generality that the numeric 
eigenproblem (2) has a reasonably accurate solution (if necessary, we may slightly per- 
turb M. and increase M, accordingly). Using ball arithmetic, it is straightforward to 
compute the matricial ball 


B(Ne, N,) = B(Ie,0)~*B(Mce, M.) B(Te, 0). 


If our numerical algorithm is accurate, then the non diagonal entries of B(.N., Ny) tend to 
be small, whence B(N., N,) can be considered as a small perturbation of a diagonal matrix. 
If we can estimate how far eigenvalues and eigenvectors of diagonal matrices can drift away 
under small perturbations, we thus obtain a solution to the original certification problem. 

Section 2 introduces notations. In Section 3, we perform a detailed study of the eigen- 
problem for small perturbations M of diagonal matrices. We exhibit a Newton iteration for 
finding the solutions. This iteration has quadratic convergence in the absence of multiple 
eigenvalues and is also an efficient tool for doubling the precision of a solution. However, 
in the case of multiple eigenvalues, the eigenproblem is ill-posed. Indeed, by a well-known 
observation, any vector occurs as the eigenvector of a small perturbation of the 2 x 2 
identity matrix. The best we can hope for is to group eigenvectors with close eigenvalues 
together in “clusters” (see also [16]) and only require T^! MT to be block diagonal. For this 
reason, we present our Newton iteration in a sufficiently general setting which encompasses 
block matrices. We will show that the iteration still admits geometric convergence for 
sufficiently small perturbations and that the blockwise certification is still sufficient for the 
computation of rigourous error bounds for the eigenvalues. In Section 4, we will present 
explicit algorithms for clustering and the overall certification problem. 


2. NOTATIONS 


2.1. Matrix norms 


Throughout this paper, we will use the max norm for vectors and the corresponding matrix 
norm. More precisely, given a vector v € C” and an n x n matrix M € C"*", we set 


lol = maxi[uil, ..., [vn|} 
|M|| = max ||Me |. 
Ivl|-1 


For a second matrix N € C"*", we clearly have 


IM NI < MIIAN || 
IMNI < MINT. 


Explicit machine computation of the matrix norm is easy using the formula 
|M| = max(|Mi| +: t [Mis 3 <in}. (3) 
In particular, when changing certain entries of a matrix M to zero, its matrix norm || M || 


can only decrease. 


2.2. Clustering 


Assume that we are given a partition 
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Such a partition will also be called a clustering and denoted by I. Two indices i, j are said 
to belong to the same cluster if there exists a k with (4, j} C Ij and we will write i + j. 
Two entries M;,; and My j’ of a matrix M € C"*" are said to belong to the same block if 
ic j and i'~ j’. We thus regard M as a generalized block matrix, for which the rows and 
columns of the blocks are not necessarily contiguous inside M. 

A matrix M € C" *" is said to be block diagonal (relative to the clustering) if M;,;—0 
whenever i^ j. Similarly, we say that M is off block diagonal if M; ; — 0 whenever i~ j. 
For a general M € C"*?, we define its block diagonal and off block diagonal projections 
A(M) = A! (M) and A(M)=Q(M) by 


Mij if inj 


_ 0 ifing 
A(M)i.3 B { 0 otherwise 


800-1 Mi, j otherwise 


By our observation at the end of section 2.1, we have 


IA) < IMI] 
JAM) < IMI. 
For the trivial clustering I, = {k}, the matrices A(M) and Q(M) are simply the diagonal 


and off diagonal projections of M. In that case we will also write A* = A and 0* =Q. 


2.3. Diagonal matrices 
Below, we will study eigenproblems for perturbations of a given diagonal matrix 
^ 
D= "n : (5) 
An 
It follows from (3) that the matrix norm u= ||D || of a diagonal matrix D is given by 
u anexos As] 

It will also be useful to define the separation number o* = o*( D) by 

g* = min{[à\;— Azji Ej} 


More generally, given a clustering as in the previous subsection, we also define the block 
separation number a = o(D) =o! (D) by 


c = min{|A;-—A,|:1~ 7} 


This number © remains high if the clustering is chosen in such a way that the indices i, j 
of any two “close” eigenvalues À; and A; belong to the same cluster. In particular, if o > 0, 
then A; — À; implies i~ j. 


3. EIGENPROBLEMS FOR PERTURBED DIAGONAL MATRICES 


3.1. The linearized equation 


Let D be a diagonal matrix (5). Given a small perturbation 


M = D+H 
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of D, where H is an off diagonal matrix, the aim of this section is to find a small matrix 
E € C?*? for which 


M' = (1+E) 1 M(14- E) 
is block diagonal. In other words, we need to solve the equation 
Q((123- E) (D+H) (1-- E)) = 0. 

When linearizing this equation in £ and H, we obtain 

O([D, E]-- H) = 0. 
If E is strongly off diagonal, then so is |D, E], and the equation further reduces to 

[D,E] = -Q(A). 

This equation can be solved using the following lemma: 


LEMMA 1. Given a matrix A € C"*" and a diagonal matrix D with entries M, ..., An, let 
B —9(D, A) € C"*" be the strongly off diagonal matrix with 


n 0 jig 
Bj = A otherwise 
Then ||B|| € a1 ||A|| and 
[D,B| = —Q(A). (6) 


Proof. The inequality follows from (3) and the definition of c. One may check (6) using 
a straightforward computation. 


3.2. The fundamental iteration 


In view of the lemma, we now consider the iteration 


(D,H) —> (D', H?), 


where 
E = OD. H 
M = (1+E) (D 4 H) (1-- E) 
D' = (M) 
H' = Q'(M!) 


In order to study the convergence of this iteration, we introduce the quantities 


B = |D| p' = J||D'| 
c = o(D) g^ = D) 
m = |A(A)| m = IAH’ 
"n, = |] m = |A(H)| 
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LEMMA 2. For ó € (0,1], assume that 


M + 72 


Then ||D' — D|| X ór and 


IN IN. VW JA 


Proof. We have 


where 
R = E(1+E) ! (D- H)(0- E)  E(D--H) E* (H, E]. 


Setting ¢=||E||<o-!m<ad< i the remainder R is bounded by 


1 
IRI < &4— a5) n(1*c) (105) ne 2(m-* m)e 
2e? 

= Toe (1 +09) H+ 2(m + 2) € 

< (4eu+2aûu)e 

< bauo im 

< n. 

Consequently, 
|D'- D|| = |[A*(M" - D)|| 2 | A*CRO| 

< |R| €óno 


m = [AGIT ISA) = [9 * CA CH + R)I 
< |H 4 R|| € m ómno 

m = IENS IRMIS IR) 
< õn. 


The inequalities u’ € u+ ôn and o/ 2 6 — 2 ôm follow from ||D' — D|| € dn. 


3.3. Convergence of the fundamental iteration 


THEOREM 3. Assume that 
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Then the sequence 
(D, H), (D', H’), (D, H^), ... 
converges geometrically to a limit (M9), HO?) with ||D9 — M|| < m and || H| 


m+n. The matriz D) + H is block diagonal and there exists a matrix Ê with |£ || 
307l m, such that 


< 
< 


D+H) = (1+É)-1(D+H)(1+É). 


Proof. Let (D), H(?) stand for the i-th fundamental iterate of (D, H) and EY = (HV), 
D®). Denote 4/9 = DO, 6 = o(D®), nË = AHO] and nf? = [AHO]. Let us 
show by induction over à that 


ID® -D| < a-5)m 
uO < uc -3)m 
o > A+) 0 
"x m (1-5) m 
np < zm. 


This is clear for i=0. Assume that the induction hypothesis holds for a given i and let 


(i) _ : ol) 1 
Q = I< —— > 
6 4 


Since (1 — z) 7 X 2 Lt, the induction hypothesis implies 


y? « 2y 
c? 2 jc 
al) > ia. 


Applying Lemma 2 for (D, HG) and ó— i we thus find 


IDD pj < [p — D] - p*9 - pj 
< (1-3) m ge E — uum 
wD < uO E20 ur -uu)m 
a) > of) Sm 2lü4l gat) 025 a 9 
mM? < nf tin? <m+ (1) m 
P" < Herr 


This completes the induction. 


Applying the induction to the sequence starting at D), we have for every j 2 0, 


aes ; 1 j 1 1 
IDOI - DO < i ea 


8 EFFICIENT CERTIFICATION OF NUMERIC SOLUTIONS TO EIGENPROBLEMS 


This shows that D? is a Cauchy sequence that tends to a limit D) with DS) — D|| € s. 
From this inequality, we also deduce that pD) — DJ < — 
geometrically to D). 


72, SO pe converges 


Moreover, for each i, we have e? = | E? || z o n = at 1». Hence, the matrix 


B=(14+ EO) (1+ EO) (14 EQO)...—1 
is well defined, and 


log(1 + || E ||) log (1 + €) + log(1 +e) + log(1 4- €) +- 


< 
< 20711. 
We deduce that 
IE « ee m-1«3o7 n, 
since g7! n2 X z 
We claim that M = D + H® converges geometrically to 
M= (14 E)! MO (14 Ê). 
For any matrix M, E € C"*" with | E|| <e « 1, we have 


(1+ £)-'M (14- E) — M || 


|IME— E(1-- E)! M (14 E)|| 


< IMI] (e+e (0 e) (1 E) HD) 
< e[MI (1 - (1-5) (1—8)73) 
2€ 
= MI (7) 


Let EO = (14+ E®) (14+ EGY) (14 EG?) ... — 1. By the same arguments as above, we 
have ;:= | EO < 307! nf) =- 07! m. Since M = (1+ £9))-1 MO (14 £9), the 
inequality (7) implies 
2 êi 
1-é; 
2 êi 


IMG) MOI < (DO || + | gr up 


= 


i (uir m nf?) 
3o m 
9i 


< 
= 1-é; 


(u + m + 72). 


This shows that M © converges geometrically to M (^). We deduce that the sequence 
H® = M) — po also converges geometrically to a limit H^) with || gr 9| < m+ me. 
Since lim;_, n? — 0, we finally observe that M (œ) = pœ) + H (99) is block diagonal. O 


THEOREM 4. Assume Ij — (k) for all k. Then, under the assumptions of Theorem 3, the 
sequence (D, H), (D', H^), (D, H"),... converges quadratically to (D(®), 0). 

Proof. The extra assumption implies that ne — 0 for all i. Let us show by induction over 2 
that we now have 
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This is clear for à = 0. Assume that the result holds for a given i. Then we may apply 
Lemma 2 to (DU), H(?) for 6=2-?'+1, and obtain 


+1 
n Pe a a 


Since ||D6*9 — DO < nC, this establishes the quadratic convergence. 


4. ALGORITHMS 


4.1. Clustering 


Let M = D + H be the perturbation of a diagonal matrix (5) as in the previous section. 
In order to apply theorem 3, we first have to find a suitable clustering (4). For a given 
threshold separation 9, we will simply take the finest clustering (i.e. for which p is maximal) 
with the property that |A; — Aj| < ô= i~ j. This clustering can be computed using the 
algorithm Cluster below. 


Algorithm Cluster 
Input: eigenvalues A1,..., An € B and 6 € IF? 
Output: the finest clustering (4) with |A; - Aj] «0 i j 


e Let G be the graph with vertices 1, ..., n and such that i and j are 
connected if and only if |A; - Aj| <ô. 


e Let H be the transitive closure of G. 
e Let M,..., Hp the connected components of H. 


e Let J; be the set of vertices of Hj, for each k. 


4.2. Certification in the case of perturbed diagonal matrices 


In order to apply theorem 3, it now remains to find a suitable threshold ó for which the 
conditions of the theorem hold. Starting with ô = 0, we will simply increase ô to o (D) 
whenever the conditions are not satisfied. This will force the number p of clusters to 
decrease by at least one at every iteration, whence the algorithm terminates. Notice that 
the workload of one iteration is O(n?), so the total running time remains bounded by O(n?). 


Algorithm DiagonalCertify 

Input: a diagonal ball matrix D € IB"*? with entries A1,..., À, and 
an off diagonal ball matrix H € B"*” 

Output: a clustering J and ê € F such that, for any M € D and H € H, the 
conditions of theorem 3 hold and ||E|| < ê 


ô:=0 

Repeat 
Compute the clustering J for A1, ..., An and ô using Cluster 
Let u:=||D||, o :— e (D), m :— A (H)|| and 9: [9 (| 


Let o :— min 4 £.,1 
6p’ 4 . 
If [m+ mal € [| and [n2] € 55, then return (J, [228]) 
Set ó:— [o] 
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4.3. 


Certification of approximate eigenvectors and eigenvalues 


Let us now return to the original problem of certifying a numerical solution to an eigen- 
problem. We will denote by 1, the n x n matrix of which all entries are one. 


Algorithm EigenvectorCertify 

Input: M —B(M., M,) € B?*" and T; c Fiij"*” 
such that p MT. is approximately diagonal 

Output: a clustering J and T= B(T., T+) € IB^*" such that for any M € M, 
there exists a T € T for which T^! MT is block diagonal 


Compute D := B(T.,0)-! M B(T}, T,) 

Let (1,&) :— DiagonalCertify( A* (ID), Q*(D)) 
Let E:— B(T.,0) B(0,&) In 

Let (Doy = LESS | for all i, j 

Return (I, B(Te, T;.)) 


Obviously, any eigenvalue A € C of a matrix M € C”*” satisfies |A| < || M |[. We may thus 
use the following modification of EigenvectorCertify in order to compute enclosures for the 
eigenvalues of M. 


Algorithm EigenvalueCertify 

Input: M —B(M,, M,) € B^*" and T; c F[i|^ *? 
such that pt MT. is approximately diagonal 

Output: ball enclosures A1,..., An € B for the eigenvalues of M, 
with the appropriate multiplicities in cases of overlapping 


Compute D := B(T.,0)-! M B(T}, T,) 
Let (1,&) :— DiagonalCertify( A* (ID), O*(.D)) 
Let m= [A (0*(D))| and 2:= [9 (9*(D))I 
For each k € {1,..., p) do 
If Ik= (i) for some i, then let A;:— B((De)i,i, [m2|) 
Otherwise 
Let c be the barycenter of the D; ; with 4 € Ik 
Let r be the maximum of |D; ; — c| for i € Ik 
Let A;:— c + B(0, [r + m - 2 ]) for all i € I, 
Return (Aj, ..., An) 


5. POSSIBLE EXTENSIONS 


Let M € C"*" be a matrix with a (numerically) multiple eigenvalue À. We have already 
stressed that it is generally impossible to provide non trivial certifications for the corre- 
sponding eigenvectors. Nevertheless, two observations should be made: 


If the eigenspace E corresponding to À has dimension 1, then small perturbations 
of the matrix M only induce small perturbations of À and E. 


Let Fy denote the full invariant subspace associated to the eigenvalue A (or all 
eigenvalues in the cluster of A). Then small perturbations of M only induce small 
perturbations of À and F). 


More precisely, in these two cases, we may search for ball enclosures for orthonormal bases 
of the vector spaces E resp. Fy, which do not contain the zero vector. 
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When considering the numeric solution (1) of the eigenproblem for M, the column 
vectors which generate P are usually far from being orthogonal. Orthonormalization can 
only be done at the expense of making T~! M T only upper triangular. Moreover, the 
orthogonalization implies a big loss of accuracy, which requires the application of a correc- 
tion method for restoring the accuracy. It seems that the fundamental Newton iteration 
from Section 3.2 can actually be used as a correction method. For instance, for small 
perturbations of the matrix 


à 1 0 0 
0 A10 0 

2 0 0 À 1 p 
0 0 0 As 


it can be shown that the fundamental iteration still converges. However, for more general 
block diagonal matrices with triangular blocks, the details are quite technical and yet to 
be worked out. 

Yet another direction for future investigations concerns the quadratic convergence. As 
a refinement of Lemma 1, we might replace D by a block diagonal matrix with entries 


A4, ..., Ap. Instead of taking B; ; = M 
jj ^ 


BijMj-MBij = Mi. 


we then have to solve equations of the form 


If the A; are sufficiently close to À; Id, it might then be possible to adapt the fundamental 
iteration accordingly so as to achieve quadratic convergence for the strongly off diag- 
onal part. 
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